function [optima, path, budget] = LocalSearch_BoxConstraint(obj,bounds,starts,step)
%LocalSearch_BoxConstraint: the local search algorithm for box constrainted
%   continueous minimization problems
%neighboring points are defined as the points where one of the dimension moves a step
% INPUTS:
%   obj: the objective function handle, f = obj(x)
%   bounds: 2 -by- d matrix, [lb;ub]
%   starts: n -by- (1+d) matrix, the start points, [objective value, x1, x2 ... xd]
%   step: the value of one dimension of the neighbor points changes
% OUTPUTS:
%   optima: n -by- (1+d) matrix, the final optimal solutions, [objective value, x1, x2 ... xd]
%   path: n -by- 1 cells, the path of the current best point
%   budget: n -by- 1 vector, the total budget used from each start
%       points (exclude start points)
[n,d] = size(starts); % the number of start points
d = d-1; % the dimension of the decision vector
optima = zeros(n,d+1);
path = cell(n,1);
budget = zeros(n,1);
for i = 1:n
    currentbest = starts(i,:);
    path{i} = currentbest;
    while 1
        % generate the neighboring points
        neighbor0 = zeros(2*d,d);
        neighbor0(1:(2*d+2):(2*d*d)) = 1;
        neighbor0(2:(2*d+2):(2*d*d)) = -1;
        NumNeighbor = size(neighbor0,1);
        neighbor0 = neighbor0.*step;
        neighbor = repmat(currentbest(2:end),[NumNeighbor,1]) + neighbor0;
        % remove points beyoud the bounds
        belongs = prod(neighbor>=repmat(bounds(1,:),[NumNeighbor,1]),2)...
            .*prod(neighbor<=repmat(bounds(2,:),[NumNeighbor,1]),2);
        neighbor(belongs==0,:) = [];
        NumNeighbor = sum(belongs);
        % objective values
        fvalue = obj(neighbor);
        budget(i) = budget(i)+NumNeighbor;
        % find the best neighboring point
        [~,best_temp] = min([currentbest(1);fvalue]);
        if best_temp==1
            break
        else
            currentbest = [fvalue(best_temp-1),neighbor(best_temp-1,:)];
            path{i} = [path{i};currentbest];
        end
    end
    optima(i,:) = currentbest;
end
end

